C This can be compiled into an xCOD resource and inserted into NumberCrunch.
C It needs an sCOD resource with the same ID, which will contain the Ascii characters:
C
C     xCOD xPROPOGATE(dT,zN:extended ;  X,Y: array ; program VXY(n:extended; xy:array) ), propogates a classical trajectory 
C
C

	SUBROUTINE PROPOGATE(dT, zN, X, Y, VXY)
	implicit EXTENDED (a-h,o-z), INTEGER (i-n)
	dimension X(*), Y(*)
	external VXY

C PROPOGATE will send a particle forward along a 2 dimensional
C classical trajectory. The gradient of the potential energy
C (which is the negative of the force) is calculated by 
C the external subroutine or user program VXY, 
C which is called via
C    Varray[1] = x
C    Varray[2] = Y
C    CALL VXY(2.0,Varray)
C  Varray has (x,y) as input to VXY, and returns (dV/dX, dV/dY) as output.
C 
C Input to xPROPOAGATE:
C			dT is the time step,
C			zN is the size of the X and Y arrays,
C			X(1), X(2), Y(1), and Y(2) are set to start the trajectory,
C			VXY is an external subroutine or user program to calculate the gradient.
C Output:
C           X(3 to zN) and Yarr(3 to zN) have been calculated.
C
C Note that 
C		- X and Y MUST point to zN component arrays (or the system crashes)
C       - (zN-3) time steps forward are calculated.
C
C The algorithm used is just:
C	X(n+1) = 2*X(n) - X(n-1) - dT^2 * Vx( X(n), Y(n) )
C	Y(n+1) = 2*Y(n) - Y(n-1) - dT^2 * Vy( X(n), Y(n) )
C
C which you get from F=ma with m=1 , F=[-Vx,-Vy] , and 
C    aX(n) = [X(n+1) - 2X(n) + X(n-1)]/dT^2, 
C which is the simplest possible discretization.
C
C						-Jim Mahoney,  7/11/89

C
C  
	extended Varray(2)
	extended Two
	
	Two = 2.0
	
	Nmin1 = zN-1
	dT2 = dT*dT
	
	do i=2,Nmin1
	
		Varray(1) = X(i)
		Varray(2) = Y(i)
		call VXY(Two, Varray)
	
		X(i+1) = 2.0*X(i) - X(i-1) - dT2 * Varray(1)
		Y(i+1) = 2.0*Y(i) - Y(i-1) - dT2 * Varray(2)		
	end do

	return
	end

C That's itÉ